#ifndef PW_PBL_FASTDIVISOR_HPP
#define PW_PBL_FASTDIVISOR_HPP

#include "IntFor.hpp"
#include "pw_pbl/Math.hpp"
#include "pw_pbl/IntFor.hpp"
#include <climits>
#include <cstdint>
#include <type_traits>
#include <utility>

namespace pw_pbl
{
	template <typename T>
	struct FastDivisor;

	template <std::signed_integral T>
	requires(sizeof(T) < sizeof(std::intmax_t)) struct FastDivisor<T>
	{
		using UT = std::make_unsigned_t<T>;
		static constexpr auto EBitCnt = sizeof(T) * CHAR_BIT;
		static constexpr auto EMaxBitPos = EBitCnt - 1;
		static constexpr UT ETwo31 = UT(1) << EMaxBitPos;

		T Magic;
		T value;
		int8_t shift;
		int8_t sign;

		explicit constexpr operator T() const { return value; }

		// divisor != 0
		// Hacker's Delight, Second Edition, Chapter 10, Integer Division By Constants
		constexpr FastDivisor(T divisor = 0) : value(divisor)
		{
			if (std::is_constant_evaluated() && divisor == 0)
				Math::_throw("divisor must not 0");
			if (value == 1)
			{
				Magic = 0;
				shift = -1;
				sign = 1;
				return;
			}
			else if (value == -1)
			{
				Magic = 0;
				shift = -1;
				sign = -1;
				return;
			}

			UT ad, anc, delta, q1, r1, q2, r2, t;
			ad = (value == 0) ? 1 : (value < 0 ? -value : value);
			t = ETwo31 + ((UT)value >> EMaxBitPos);
			anc = t - 1 - t % ad;
			q1 = ETwo31 / anc;
			r1 = ETwo31 - q1 * anc;
			q2 = ETwo31 / ad;
			r2 = ETwo31 - q2 * ad;
			int_fast8_t p = EMaxBitPos;
			do
			{
				++p;
				q1 = 2 * q1;
				r1 = 2 * r1;
				if (r1 >= anc)
				{
					++q1;
					r1 -= anc;
				}
				q2 = 2 * q2;
				r2 = 2 * r2;
				if (r2 >= ad)
				{
					++q2;
					r2 -= ad;
				}
				delta = ad - r2;
			} while (q1 < delta || (q1 == delta && r1 == 0));
			Magic = q2 + 1;
			if (value < 0)
				Magic = -Magic;
			shift = p - EBitCnt;

			if ((value > 0) && (Magic < 0))
				sign = 1;
			else if ((value < 0) && (Magic > 0))
				sign = -1;
			else
				sign = 0;
		}

	public:
		static constexpr T calc(T n, const FastDivisor &divisor)
		{
			T q = UIntFor::Size<sizeof(T) * 2>(divisor.Magic * (SIntFor::Size<sizeof(T) * 2>)n) >> EBitCnt;
			q += n * divisor.sign;
			if (divisor.shift >= 0)
			{
				q >>= divisor.shift; // we rely on this to be implemented as arithmetic shift
				q += (UT)q >> EMaxBitPos;
			}
			return q;
		}
	};

	/// floor(n/divisor) = floor( ceil(2**P/divisor) * n / 2**P ) = floor( floor(((ceil(2**P/divisor) - 2**W) * n) / 2**W) + n / 2**P )
	/// @note W = bitwidth, P = W + log2ceil(divisor)
	/// @note Integer Division Using Reciprocals by Robert Alverson
	template <std::unsigned_integral T, T DIVISOR>
	requires(sizeof(T) < sizeof(std::uintmax_t)) struct FastDivisor<std::integral_constant<T, DIVISOR>>
	{
		static constexpr auto EBitWidth = Math::BitCntOf<T>;
		static constexpr auto ELog2DivisorCeil = Math::Log2Ceil(DIVISOR);
		static constexpr auto EShiftNum = EBitWidth + ELog2DivisorCeil;
		using TCalcUInt = UIntFor::Size<sizeof(T) * 2>;
		static constexpr T EMagicNum = {(DIVISOR - 1 + (TCalcUInt(1) << EShiftNum)) / DIVISOR - (TCalcUInt(1) << EBitWidth)};
		static constexpr T calc(T n, const FastDivisor &) noexcept
		{
			return ((TCalcUInt(EMagicNum) * n >> EBitWidth) + n) >> ELog2DivisorCeil;
		}
	};

	/// floor(n/divisor) = floor( ceil(2**P/divisor) * n / 2**P )
	/// @note W = bitwidth, P = W + log2ceil(divisor)
	/// @note Integer Division Using Reciprocals by Robert Alverson
	/// @note 当被除数和除数位宽小于类型位宽时的优化计算
	template <std::unsigned_integral T, T DIVISOR, T BITWIDTH>
	requires(sizeof(T) < sizeof(std::uintmax_t) && BITWIDTH < Math::BitCntOf<T>) struct FastDivisor<std::integer_sequence<T, DIVISOR, BITWIDTH>>
	{
		static constexpr auto EBitWidth = BITWIDTH;
		static constexpr auto ELog2DivisorCeil = Math::Log2Ceil(DIVISOR);
		static_assert(ELog2DivisorCeil <= EBitWidth);
		static constexpr auto EShiftNum = EBitWidth + ELog2DivisorCeil;
		using TCalcUInt = UIntFor::BitCnt<EBitWidth * 2>;
		static constexpr T EMagicNum = {(DIVISOR - 1 + (TCalcUInt(1) << EShiftNum)) / DIVISOR};
		static constexpr T calc(T n, const FastDivisor &) noexcept
		{
			return TCalcUInt(EMagicNum) * n >> EShiftNum;
		}
	};

	/// floor(n/divisor) = floor( ceil(2**P/divisor) * n / 2**P )
	/// @note W = bitwidth, P = W + log2ceil(divisor)
	/// @note Integer Division Using Reciprocals by Robert Alverson
	/// @note 当被除数和除数位宽小于类型位宽时的优化计算
	template <uintmax_t DIVISOR, uintmax_t BITWIDTH>
	requires(BITWIDTH < Math::BitCntOf<uintmax_t>) struct FastDivisor<std::integer_sequence<uintmax_t, DIVISOR, BITWIDTH>>
	{
		static constexpr auto EBitWidth = BITWIDTH;
		static constexpr auto ELog2DivisorCeil = Math::Log2Ceil(DIVISOR);
		static_assert(ELog2DivisorCeil <= EBitWidth);
		static constexpr auto EShiftNum = EBitWidth + ELog2DivisorCeil;
		// static_assert(EShiftNum < Math::BitCntOf<uintmax_t>);
		static constexpr uintmax_t EMagicNum = {(DIVISOR - 1 + (uintmax_t(1) << EShiftNum)) / DIVISOR};
		static constexpr uintmax_t calc(uintmax_t n, const FastDivisor &) noexcept
		{
			return EMagicNum * n >> EShiftNum;
		}
	};

	template <std::signed_integral S, std::signed_integral T>
	requires(sizeof(S) <= sizeof(T)) constexpr T operator/(S n, const FastDivisor<T> &divisor) noexcept { return FastDivisor<T>::calc(n, divisor); }

	template <std::signed_integral S, std::signed_integral T>
	requires(sizeof(S) <= sizeof(T)) constexpr T operator%(S n, const FastDivisor<T> &divisor) noexcept { return n - FastDivisor<T>::calc(n, divisor) * divisor; }
} // namespace pw_pbl

#endif // PW_PBL_FASTDIVISOR_HPP
